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Abstract The aim of the present paper is to develop a strategy for solving reliability-based design optimization 
(RBDO) problems that remains applicable when the performance models are expensive to evaluate. Starting with 
the premise that simulation-based approaches are not affordable for such problems, and that the most-probable- 
failure-point-based approaches do not permit to quantify the error on the estimation of the failure probability, an 
approach based on both metamodels and advanced simulation techniques is explored. The kriging metamodeling 
technique is chosen in order to surrogate the performance functions because it allows one to genuinely quantify the 
surrogate error. The surrogate error onto the limit-state surfaces is propagated to the failure probabilities estimates 
in order to provide an empirical error measure. This error is then sequentially reduced by means of a population- 
based adaptive refinement technique until the kriging surrogates are accurate enough for reliability analysis. This 
original refinement strategy makes it possible to add several observations in the design of experiments at the same 
time. Reliability and reliability sensitivity analyses are performed by means of the subset simulation technique for 
the sake of numerical efficiency. The adaptive surrogate-based strategy for reliability estimation is finally involved 
into a classical gradient-based optimization algorithm in order to solve the RBDO problem. The kriging surrogates 
are built in a so-called augmented reliability space thus making them reusable from one nested RBDO iteration to 
the other. The strategy is compared to other approaches available in the literature on three academic examples in 
the field of structural mechanics. 

Keywords reliability-based design optimization (RBDO) • kriging • surrogate modeling • subset simulation • 
adaptive refinement 



1 Introduction 

In structural mechanics, design optimization is the decision-making process that aims at finding the best set of 
design variables which minimizes some cost model while satisfying some performance requirements. Due to the 
inconsistency between these two objectives, the optimal solutions often lie on the boundaries of the admissible 
space. Thus, these solutions are rather sensitive to uncertainty either in the parameters (aleatory) or in the models 
themselves (epistemic) . Reliability-based design optimization (RBDO) is a concept that accounts for uncertainty all 
along the optimization process. Basically, the deterministic performance model is wrapped into a more realistic 
probabilistic constraint which is referred to as the failure probability. Despite its attractive formulation, the appli- 
cation field of RBDO is still limited to academic examples. This is mostly due to the fact that it is either based 
on simplifying assumptions that might not hold in practice; or in contrast, it requires computationally intensive 
stochastic simulations that are not affordable for real industrial problems. The present work attempts to propose 
an efficient strategy that would in fine bring the RBDO application field to more sophisticated examples, closer to 
real engineering cases. In other words, the challenge is (i) to provide an optimal safe design within a few hundred 
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evaluations of the performance models and (ii) to be able to quantify and minimize the errors induced by the 
various assumptions that are made along the development of the resolution strategy. 

The remaining part of this introduction is devoted to the formulation of the RBDO problem the authors attempt 
to solve. A short literature review is also provided as an argument for the presently proposed surrogate-based 
resolution strategy. Section 2 introduces the kriging surrogate model. A specific emphasis is put on the epistemic 
nature of the prediction error that is then used in Section 3 in order to quantify and minimize the surrogate error. 
Section 4 involves the adaptive refinement strategy of the kriging surrogate into a nested reliability-based design 
optimization loop. The convergence of the approach is finally heuristically demonstrated in Section 5 through a few 
academic examples from the RBDO literature. 



1 . 1 Problem formulation 

Given a parametric model for the random vector X describing the environment of the system to be designed, the 
most basic formulation for the RBDO problem reads as follows: 



In this formulation, c is the objective function to be minimized with respect to the design variables e 9 e , 
while satisfying to n c deterministic soft constraints {/ ; , i = 1, ... ,n c } bounding the so-called admissible design 
space defined by the analyst. Note that in most applications these soft constraints consist in simple analytical 
functions that prevent the optimization algorithm from exploring regions of the design space that have no physical 
meaning (e.g. negative or infinite dimensions), so that these constraints are inexpensive to evaluate. A deterministic 
design optimization (DDO) problem would simply require additional performance functions {g h I = 1, ... ,n p ] 
describing system failure with respect to the specific code of practice. As opposed to the previous soft constraints, 
these functions often involve the output of an expensive-to-evaluate black-box function M - e.g. a finite element 
model. RBDO differs from DDO in the sense that these constraints are wrapped into n p probabilistic constraints 
{P (gi < 0) < Pf t , I = 1, . . . , n p }. P° ( is the minimum safety requirement expressed here in the form of 

an acceptable probability of failure which may be different for each performance function g l . Such probabilities of 
failure are conveniently defined in terms of the following multidimensional integrals: 



One should notice that, in the present formulation, the design vector 9 is a set of hyperparameters defining 
the random vector X. In other words, in this work, design variables are exclusively considered as hyperparameters 
in the joint probability density function f x of the random vector X because it will later simplify the computation 
of reliability sensitivity - i.e. the gradients of the failure probability. There is however no loss of generality since 
deterministic design variables might possibly be considered as artificially random (either normal or uniform) with 
small variance - i.e. sufficiently close to zero. 

One could possibly argue that this formulation lacks full probabilistic consideration because the cost function 
is defined in a deterministic manner as it only depends on the hyperparameters of the random vector X. A more 
realistic formulation should eventually account for the randomness of the cost function possibly induced by the 
one in X. However, the present formulation is extensively used in the RBDO literature for simplicity. Note however 
that thanks to the rather low complexity of usual cost models (analytical functions), an accurate simulation-based 
estimation of a mean cost, say c (9) = E [c (X (0))] would not require a large computational effort. 



1.2 Short literature review 

The most straightforward approach to solve the RBDO problem in Eq. (1) consists in nesting a reliability analysis 
within a nonlinear constrained optimization loop. Such methods are referred to as double-loop or nested approaches. 
Despite their conceptual simplicity these approaches are often argued to lack efficiency since they require too many 
evaluations of the performance functions g t that might involve the output of a time-consuming computer code. 
However for a broad range of applications where the performance functions are linear or weakly nonlinear - i.e. 
when the first order reliability method (FORM) is applicable, the nested approach is able to give results within 
a reasonable number of evaluations of the performance functions (e.g. Enevoldsen and Sorensen (1994) achieve 
convergence within a few thousands evaluations) . 




(1) 




(2) 
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The formulation in Eq. (1) is known as the reliability index approach (RIA). Note that RIA often refers to 
nested RBDO algorithm based on FORM despite the fact that the concept can easily be extended to simulation- 
based reliability methods. Some authors, starting with Tu et al. (1999), proposed an alternative formulation to the 
RBDO problem known as the performance measure approach (PMA). The probabilistic constraints are transformed 
into quantile constraints that are approximated using the first-order reliability theory. According to Youn and Choi 
(2004), PMA would be more stable and efficient than RIA because the inner reliability algorithms used to solve the 
first-order quantile approximation (the so-called Mean Value algorithms) are argued to be much more efficient than 
their probability approximation counterparts (.e.g. the Hasofer-Lind-Rackwitz-Fiessler algorithm). 

The nested approach becomes intractable in the case of more complex performance functions for which the 
nested reliability analysis should resort to simulation-based methods. For such cases however, Royset and Polak 
(2004a,b) proposed the so-called sample average approximation method which consists in gradually refining the 
simulation-based reliability analysis as the optimization algorithm converges towards an optimal design. To do so, 
they propose an empirical stepwise refinement criteria to define whether the number of simulations should be 
raised or not. Note also that the use of variance reduction techniques can significantly reduce the overall number 
of simulation runs such as demonstrated in Jensen et al. (2009) where the authors used subset simulation. How- 
ever, simulation-based approaches still require too many performance functions evaluations to make the approach 
applicable to real engineering cases even though their application field is growing with the increasing availability 
of high performance computational (HPC) resources (e.g. interconnected clusters of PCs). 

An alternative to double-loop approach consists in decoupling the optimization loop from the reliability analyses 
so that both can be sequentially performed in an independent manner. Such approaches are referred to as sequential 
approaches or decoupled approaches (Royset et al., 2001; Du and Chen, 2004; Aoues and Chateauneuf, 2010). An 
advantage of these approaches is that they do not require any reliability sensitivity analysis since the optimization is 
performed directly onto the performance function. However, the decoupling often relies on the most probable failure 
point (MPFP) assumptions and thus suffers from the possible non-unicity of this point and the strong nonlinearities 
in the performance functions. 

Single-loop approaches (Kuschel and Rackwitz, 1997; Kirjner-Neto et al., 1998; Kharmanda et al., 2002; Shan 
and Wang, 2008) attempt to fully reformulate the original RBDO problem into an equivalent DDO problem that 
allows a simple and efficient resolution by means of classical optimization algorithms. The approach is mainly based 
on concepts that are closely related to the notion of partial safety factors. It is certainly the most computationally 
efficient approach as soon as the assumptions under which the probabilistic-deterministic equivalence is built hold. 
Once again, most of them are based upon the assumption that the MPFP exists and that it is unique. 

Stochastic subset optimization (SSO) is a simulation-based approach recently proposed by Taflanidis and Beck 
(2009a) that consists in finding the region of the admissible design space where the failure probability density 
function is minimal. It is based on conditional simulations in a so-called augmented reliability space where the design 
variables are artificially considered as random with uniform distribution. The range of uncertainty in the design 
variables is reduced along with the identification of low failure probability density regions. The overall concept is 
closely related with the subset simulation reliability method proposed by Au and Beck (2001), and further explored 
by Au (2005) for reliability-based design sensitivity analysis. In the end, the algorithm provides a set of parameters 
that is likely to contain the optimal solution which can possibly be found by means of a more refined stochastic 
search algorithm. However, the problem that SSO attempts to solve is not a full RBDO problem in the sense that it 
is designed to minimize the failure probability whereas the purpose of RBDO is to minimize a cost function under 
some maximal failure probability constraint. 

Finally, the approach that is investigated here is the surrogate-based (or response-surface-based) approach. Start- 
ing with the premise that the performance function evaluation might involve a time-consuming computational task, 
the approach consists in replacing this function by a surrogate that is much faster to evaluate. This surrogate is usu- 
ally built and possibly refined from a few evaluations of the real performance function. Various surrogates have 
been used in the structural reliability literature out of which: polynomial response surfaces (Faravelli, 1989; Bucher 
and Bourgund, 1990; Kim and Na, 1997; Das and Zheng, 2000), polynomial chaos expansions (Sudret and Der 
Kiureghian, 2002; Berveiller et al., 2006; Blatman and Sudret, 2008, 2010), support vector machines (Hurtado, 
2004; Deheeger, 2008; Bourinet et al, 2010), neural networks (Papadrakakis and Lagaros, 2002) and kriging (Bi- 
chon et al., 2008). Such surrogates are argued here to be more flexible than the classical Taylor-expansion-based 
methods (i.e. approaches using the first- or second-order reliability theory) that are commonly used in the RBDO 
literature. An additional interesting fact about this approach is that some of these surrogates, namely kriging and 
support vector machines, allow one to quantify an empirical surrogate error measure that can be propagated to the 
final quantity of interest: the failure probability. 
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2 The kriging surrogate 

In computer experiments, surrogate modeling addresses the problem of replacing a time-consuming mathematical 
model M called a simulator by an emulator that is much faster to evaluate. The emulator has the same input space 
@ x c k" and output space 9 y c| than the simulator. For the sake of simplicity only scalar-output simulators are 
considered here. It is important to note here that both x and y are deterministic in this context. In the present 
RBDO application, the time-consuming models are the performance functions {g lt I = 1, . . . , n p } in Eq. (1). 

Kriging (Santner et al., 2003) is one particular emulator that is able to give a probabilistic response Y(x) 
whose variance (spread) depends on the quantity of available knowledge. In other words, the uncertainty in this 
prediction is purely epistemic and due to a lack of knowledge at specific input x . For a more detailed discussion 
on the distinction that should be made between aleatory and epistemic sources of uncertainty, the reader is referred 
to Der Kiureghian and Ditlevsen (2009). This interesting property of the kriging surrogate will be used in the next 
sections to make the surrogate-based RBDO approach adaptive. 

In essence, kriging starts with the assumption that the output y = M [x ) is a sample path of a Gaussian stochas- 
tic process Y whose mean and autocovariance functions are unknown and to be determined from the values of the 
model response 'W = {y t = M , i = 1, ... ,m] evaluated onto an experimental design X = {x 1 , ... ,x m }. 
More specifically, the so-called universal kriging (UK) model assumes the Gaussian process has a stationary autoco- 
variance function and is centered with respect to a non-stationary trend in the form of a linear regression model so 
that it reads: 

p 

Y(x) = Y J PifiW + ZM=f{x) J p+Z{x). (3) 

i=l 

The first term is the regression part that requires the choice of a functional set /ei? 2 ($) x , M.) . The second term is 
a stationary Gaussian process with zero mean, constant variance a\ and stationary autocorrelation function R so 
that its autocovariance function reads as follows: 

C(x,x')=cr*R(\x-x'\,i). (4) 

The parameters I, /3 and shall be determined according to a given class of symmetric positive definite autocor- 
relation functions and a given set of regression functions. The most widely used class of autocorrelation functions 
is the anisotropic generalized exponential model: 

R(|x-x'|,l)=exp^- l**~* fc l j, 1<5<2. (5) 

Generally speaking, the choice of the autocorrelation model should be made consistently with the known proper- 
ties of M such as derivability or periodicity (see e.g. Rasmussen and Williams, 2006, ch. 4). Common regression 
models involve zero, first and second-order polynomials of {x k , k = 1, . . ., n}. Note however that as for ordinary 
least-squares regression, the size p of the functional basis together with the dimension of the input vector n dra- 
matically increases the number m of required observations. When attempting to model the output of a high-fidelity 
computer model, a good choice for the regression model may be an equivalent low-fidelity analytical model built 
on simplifying assumptions. 

The construction of a kriging model consists in a two-stage framework that is further detailed in Subsections 
2.1 and 2.3. 



2.1 The conditional distribution of the prediction 

The first stage aims at looking for the distribution of some prediction Y(x) = Y(x) \Y of the Gaussian process (GP) 
at any location x e Qi x conditional on the vector of observations Y = (Y [x^) , ... , Y (x m )} T . Since Y is assumed 
to be a Gaussian process, the vector of observations Y is distributed as follows: 

Y~Jf(Fp, CT^R) (6) 

where we have introduced the matrices F and R that are defined with respect to the assumed statistics of the 
Gaussian process. Their terms read as follows: 

F U = fj O;) > i = l, ■■■ ,m, j = 1, ... ,p (7) 
R U = R {\ x i- X j\' i ) ' i = l, ■■■ ,m, j = 1, ... ,m. (8) 
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It is also known that the augmented random vector (y t , Y (x)) T is jointly Gaussian by definition: 



Y 1 .,/ F 



^ tr^T P, o 



R r O) 
r(x) T 1 



(9) 



where we have introduced the vector of cross-correlations r(x) between the observations and the prediction: 

r i (x)=R(]x-x i \,e), i = l,...,m. (10) 

It is well known (see e.g. Severini, 2005, Chap. 8) that the conditional distribution of Y(x) = Y(x) | Y is also 
Gaussian with mean: 

H f (x)=f(x) T p + r(x) J R- 1 (Y-¥p) (11) 

and variance: 

<4 00 = o\ (l - r 00 T R- 1 r(x) + u(x) T (F T R" 1 F)" 1 u00) (12) 



where we have introduced: 



P — (F T R _1 F) 1 F T R 1 Y (13) 
u(x) = ¥ T R- 1 r(x)-f(x). (14) 



2.2 Properties of the kriging surrogate 



The kriging surrogate is an exact interpolator. Indeed, noting that the i-th column of the matrix R is equal to the 
vector r (x,) = e,, the following relation holds: 

Re,- = r 0,) (15) 

where e i is the i-th basis vector of K m which has all-zero components except its i-th component equal to one. Using 
the symmetry property of R and plugging the latter relation in the predictor's mean in Eq. (11) and in its variance 
in Eq. (12) finally allows to prove the following statement: 

i = 1 '-' m ' (16) 

which means that the kriging surrogate interpolates the observations with variance equal to zero (deterministic 
prediction) . This property holds for any regression and regular autocorrelation models, whatever the values of the 
parameters I, p and a y . 

2.3 Joint maximum likelihood estimation of the Gaussian process parameters 

In the first stage it was assumed that the autocovariance function u Y R{»,t) and the regression model / («) T p are 
known. In computer experiments or in geostatistics it is never the case so that the models need to be chosen and their 
parameters must be determined using common statistical inference techniques such as the variographic analysis 
(VA), maximum likelihood estimation (MLE) or Bayesian estimation (BE). 

Historically, the kriging methodology was introduced by geostatisticians in order to predict a map of soil proper- 
ties from in situ observations for mining prospection. In this field, the autocovariance structure is usually estimated 
from the data using variographic analysis (Cressie, 1993). Then, provided the empirical variogram features some 
required properties it can be turned into an autocovariance model. However, this methodology is not well-suited to 
our purpose because of its user-interactivity. 

In computer experiments, the most widely used methodology is the MLE technique. Provided a functional set 
/ e S£ 2 (S) x , R) and a stationary autocorrelation model R(», I) are chosen, one can express the likelihood of the 
data with respect to the model and maximize it with respect to the sought parameters (i, p and a Y ). One can 
show that p and a Y can be derived analytically (using the first-order optimality conditions) and solely depends 
on the autocovariance parameters £ that are solution of a numerically tractable global optimization problem - see 
e.g. Welch et al. (1992); Lophaven et al. (2002) for more details. This technique, implemented within the DACE 
toolbox by Lophaven et al. (2002), was used for the applications presented in this paper. 

Note that the BE technique consists in weighting the likelihood by a prior over the autocovariance parameters 
and finally allows to derive a full probabilistic posterior on them. Despites it has an attractive formulation as it 
allows one to quantify an additional epistemic uncertainty in the autocovariance parameters themselves, it is still 



6 



V Dubourg et al. 



difficult to involve it in a nested general-purpose procedure such as the one proposed in this paper due to the 
required prior information which might not be available. 

In the case of MLE, the variance o~| underestimates the real variance as it does not account for (i) the empirical 
choice of the autocovariance structure of the GP, and (ii) the uncertainty induced by the statistical inference of 
I . Note that BE does account for this uncertainty but it is rather difficult to propagate this additional uncertainty 
through the predictor - this is sometimes referred to as Bayesian kriging. Even though, the probabilistic response of 
the kriging model is convenient as it allows one to: 

- provide probabilistic confidence intervals in addition to the prediction; 

- define refinement criteria to build the kriging surrogate in an adaptive refinement procedure depending on the 
region of interest as introduced in Section 3. 



3 Design of experiments for the kriging surrogate 

Various refinement techniques have been proposed in the kriging-related literature, e.g. for global optimization 
(Jones et al., 1998) or for probability/quantile estimation (Oakley, 2004; Bichon et al, 2008; Lee and Jung, 2008; 
Ranjan et al., 2008; Vazquez and Beet, 2009; Picheny et al., 2010). They are all based on the genuine probabilistic 
nature of the kriging prediction. An adaptive refinement technique consists in finding the set of points that should 
be added to the design of experiments (DOE) in order to improve the prediction and reduce its associated epistemic 
uncertainty in some specific region of interest (Picheny et al., 2010). In global optimization, the region of interest is 
the vicinity of the optimum that can be sequentially achieved by means of the so-called efficient global optimization 
technique (Jones et al., 1998). In reliability and reliability-based design optimization, the region of interest in the 
space of random variables X is the vicinity of the limit-state surface g(x ) = (we drop the subscript I in this section 
for the sake of clarity) . 



3.1 Premise 

Provided an initial kriging prediction of the performance function g denoted by G ~ jV (jug, erg), the probable 
vicinity of the limit-state surface g = can be defined in terms of the following margin of uncertainty (epistemic 
uncertainty) : 

M={x : -ka 5 (x)<u 5 (x)<+ka 5 (x)} (17) 

where k might be chosen as k = <P^ (97.5%) = 1.96 meaning that a 95% confidence interval onto the prediction 
of the limit-state surface is chosen. Any point that belongs to this margin correspond to an uncertain sign of the 
prediction. The spread of this margin should be reduced in order to make the reliability estimation accurate, and 
this can be achieved by adding new points in the DOE that are located in this margin. 

The probability that a point x e @ x falls into the margin M (with respect to the epistemic uncertainty in the 
prediction) can be expressed in closed-form as follows: 

fkor (x) — up (x)\ f —kap(x) — Ur(x)\ 

P(xeM) = <i> — , 7^ -<£ - . (18) 

V °gM J V °gM J 

Note that the same criterion is used in Picheny et al. (2010) where it has been named the weighted integrated 
mean squared error (W-IMSE). Then, finding the point that maximizes this quantity on the support of the PDF of 
X will finally bring the best improvement point in the DOE. Starting with this statement, many authors in the 
kriging literature decide to use global optimization algorithms in order to find the best improvement point. For 
instance, Bichon et al. (2008) use a different criterion named the expected feasibility function, and Lee and Jung 
(2008) use the constraint boundary sampling criterion, but both use a global optimization algorithm to find the best 
improvement point. 



3.2 Principle of the proposed refinement procedure 

In this paper, we propose a different strategy that allows one to add a set of improvement points. Indeed, all the 
proposed refinement criteria (including the one in Eq. (18)) are highly multimodal, thus meaning that: 

- the global optimization problem is not easy to solve; 

- there does not really exist one best point (especially for the contour approximation problem of interest), since 
each mode is associated with a potentially interesting point to add to the DOE. 
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It may also be faster if one has the ability to perform several simulation runs of g at the same time (e.g. on some 
distributed computational platform). Note that the idea is inspired from Hurtado (2004); Deheeger and Lemaire 
(2007); Deheeger (2008); Bourinet et al. (2010) where the authors use the same idea applied to support vector 
classifiers. 

To do so, the probability to fall in the margin of epistemic uncertainty is multiplied by a weighting PDF w so 
that the following criterion: 

^(x) = P(xeM)w(x) (19) 

can itself be considered as a PDF up to an unknown but finite normalizing constant. The weighting density w can 
either be the original PDF of X, or, as it is proposed here, the uniform PDF on a sufficiently large confidence region 
of the original PDF. 

Such a confidence region might be difficult to define for any given PDF, but as it is usually done in structural 
reliability (Ditlevsen and Madsen, 1996), the original random vector X can be transformed into a probabilistically 
equivalent standard Gaussian random vector U for which the confidence region is simply an hypersphere with 
radius p . The reader is referred to Lebrun and Dutfoy (2009) for a recent discussion on such mappings u = T(x). 
In that given space, p o can be easily selected as e.g. p = 8 which corresponds to the maximal reliability index 
that can be justified numerically, and the sought uniform PDF is simply defined in terms of the following indicator 
function: 

w ^ Kl v^W u) - (20) 

The normalizing constant of this PDF is the volume of the hypersphere and it could thus be easily derived though 
it is not required to generate samples from ^ as explained hereafter. Note that another weighting PDF will be 
introduced later on (in Section 4.1) for the RBDO problem of interest. 

Given that pseudo-PDF, one can generate JV (say JV = 10 4 ) samples by means of any well suited Markov chain 
Monte Carlo (MCMC) simulation technique (Robert and Casella, 2004), such as the slice sampling technique (Neal, 
2003). The N generated candidates are expected to be highly concentrated at the modes of ^ - and thus in M. 
This population of candidates is then reduced to a smaller one that has essentially the same statistical properties 
(i.e. it uniformally covers M) by means of the K -means clustering technique (MacQueen, 1967). 

We make a final remark considering the three approximate failure subsets F _1 , F °, F +1 defined as follows: 



F' = {x : /xg(x) + ifcc7g(x)<0}, i = -l,0,+l. (21) 

Note that the extreme failure subsets are related with the margin of uncertainty through the following set relation- 
ship: M = F -1 \ F +1 . Due to the positiveness of standard deviation, the following statement holds: 

F+ 1 C F C F" 1 => P (x G F +1 ) < P (x G F °) < P (X G F" 1 ) . (22) 

Unfortunately, there is no proof that the real failure subset F satisfies the following relationship F +1 ck F -1 due 
to the empirical assumptions that were made in Section 2. However, the spread of the interval [P(X g F +1 ), P(X g 
F -1 )] is a useful measure to check if the kriging surrogate is accurate enough for reliability analysis or not, and it is 
used in this paper as a stopping criterion for the proposed refinement procedure. Note that due to the inconvenient 
order of magnitude of low probabilities, it is more meaningful to work with the generalized reliability indices that 
are defined as follows: 

p ' = —<P^ 1 ^P(X 6F')J, i = -l, 0, +1. (23) 

Then, in the proposed applications, the accuracy criterion is usually set to = 10 1 — 10~ 2 and the refinement 
stops when: 

max (j3 +1 - /3 °, p -^ 1 ) <e p . (24) 

When p 1 are estimated by means of simulation techniques (Monte-Carlo or subset simulation), the proposed 
accuracy criterion should account for the additional uncertainty induced by the lack of simulations. To do so, in the 
present paper, p^ 1 and p +1 estimates are replaced with their respective lower and upper 95% confidence bounds 
based on their associated variance of estimation. Thus, should be selected in accordance with the given number 
of simulations used for reliability estimation. 
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3.3 Implementation 

In order to summarize the proposed refinement procedure, we provide the pseudo-code in Algorithm 1. First, we 
initialize the empty DOE [3C = %, ( S = 0), the uniform refinement pseudo-PDF ^ for the first space-filling DOE 
and we select the level of confidence k in the metamodel. Then, we generate the candidate population 3? from 
the density function by means of any well-suited MCMC simulation technique (using e.g. slice sampling). 
This population is reduced to its K clusters' center using K-means clustering - K being given. The performance 
function <& is evaluated onto these K newly selected points X new and a new kriging model is built from the updated 
DOE {% , Note that the kriging model construction step involves the maximum likelihood estimation of the 
autocorrelation parameters I £ [li,£u] - The refinement pseudo-PDF V is also updated. Finally, a reliability analysis 
(using e.g. subset simulation) is performed onto the three approximate failure subsets F +1 , F °, and F -1 in order to 
compute the proposed error measure. The DOE is enriched if and while this error measure exceeds a given tolerance 
Ep = KT 1 - 1(T 2 . 



Algorithm 1 Population-based adaptive refinement strategy 
_______ _ 

2: w := x — > w(x) 

3: k :=$- 1 (97.5%), s fj := 10" 1 

4: Refine := true 

5: ¥(x 6M) :=x ~ 1 

6: while Refine do 

7: & := MCMCAlgorithm(^) 

8: 3C nevl := KMeansAlgorithm^ , K) 

9. ^new ■ S (*^new) 

io: sc -.= {s:,s: new }, <§:={<§, % evi } 

11: 6 := MaximumLikelihoodKrigingModelCr, »,/(•),«(•, I), I L , t u ) 

12: Pfr6 M3«^#(ia^)_ # (_^__) 

13: <€ :=x ^¥(x 6i) w(x) 

14: for i:= -1,0,-1-1 do 

15: F i :={x : fig(x) + i fccrg(x) < o} 

16: p 1 := ReliabilityAnalysis(F ') 

17: end for 

18: Refine := max (/3 +1 - /S °, P - P' 1 ) > e p 

19: end while 



Figure 1 illustrates the proposed adaptive refinement strategy applied to a nonlinear limit-state surface from 
Waarts (2000) . The upper subfigures show the contours of the refinement pseudo-PDF at each refinement step 
together with the candidate population generated by slice sampling and its K clusters' center obtained by if -means 
clustering - K = 10 being given. For this application, the weighting PDF w was selected as the uniform density 
in the /3 -radius hypersphere - /3 = 8. It can be observed that the refinement criterion features several modes as 
argued earlier in this Section. In the lower subfigures one can see the real limit-state surface represented as the 
dashed black curve, its kriging prediction represented as the black line and its associated margin of uncertainty M 
which is bounded below by the red line and above by the blue line. Another interpretation of these figures is that 
any point within the blue bounded shape is positive with a 95% confidence and any point inside the red bounded 
shape is negative with the same confidence level. 



4 The proposed adaptive surrogate-based RBDO strategy 

The proposed strategy to solve the RBDO problem in Eq. (1) consists in nesting the previously introduced kriging 
surrogate together with the proposed refinement strategy within a classical (but efficient) nested RBDO algorithm. 



4.1 The augmented reliability space 

In this section, we describe the space where the kriging surrogates are built. Indeed, observing that building the 
kriging surrogates from an empty DOE for each nested reliability analysis [e.g. in the space of the standard normal 
random variables) would be particularly inefficient, it is proposed to build and refine one unique global kriging 
surrogate for all the nested reliability analyses. 
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"1 "1 1 

(a) Init.: 10-point DOE (b) Step 3: 40-point DOE (<0 Ste P 9: 100-point DOE 



Fig. 1 An illustration of the proposed limit-state surface refinement technique 



Such a globality can be achieved by working in the so-called augmented reliability space such as defined in 
Taflanidis and Beck (2009b). In Kharmanda et al. (2002), the augmented reliability space is defined as the tensor 
product between the space of the standardized normal random variables and the design space: Q v x @ e , but the 
dimension of this space (n + n e ) suffers from both the number of random variables n and the number of design 
variables n e . It is also argued here that this space may cause some loss of information as the performance functions 
are not in bijection with that augmented space. In contrast, in Taflanidis and Beck (2009b) and in the present 
approach, the dimension of the augmented reliability space is kept equal to n by considering that the design vector 
9 simply augments the uncertainty in the random vector X. Indeed, the augmented random vector V = X (0) has 
a PDF h which accounts for both an instrumental uncertainty in the design choices & and the aleatory uncertainty 
in the random vector X. Under such considerations, h reads as follows: 

Hv)= I f x (x\0)n(6)d0 (25) 
J% 

where f x is the PDF of X given the parameters and n is the PDF of that can be assumed uniform on the design 
space @ e . An illustration of this augmented PDF is provided in Figure 2 in the univariate case. The augmented 
reliability space is spanned by the axis V on the left in this simple case. 

The DOE should cover uniformally a sufficiently large confidence region of this augmented PDF in order to 
make the surrogate limit-state surfaces accurate wherever they can potentially be evaluated along the optimization 
process. More precisely, they should be accurate for extreme design choices 9 (i.e. located onto the boundaries of 
the optimization space @ e ) and extreme values of the marginal random vector X (to be able to compute reliability 
indices as large as e.g. /3 = 8). A confidence region is essentially the multivariate extension of the univariate 
concept of confidence interval. Under the previous general assumptions, it is hard to give a mathematical form to 
the contour of this region. However, one may easily build an hyperrectangular region that bounds the confidence 
region of interest. 

Indeed, such an hyperrectangular region is defined as the tensor product of the confidence intervals on the 
augmented margins {V t , i = 1, . . . , n}. In order to compute the quantiles bounding these confidence intervals, one 
should additionally assume that the design parameters are exclusively involved in the definition of the margins 
- i.e. no parameters in the dependence structure (the copula) as it will never be the case in most RBDO applica- 
tions. For each margin, the lower quantiles {q^ , i=l, ... ,n] (at the probability level 1 — $ (/3 ) = # (— fti)) an d 
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Fig. 2 The augmented probability density function of a Gaussian random variate with uniform mean. 

upper quantiles {qj, i = 1, ... ,n} (at the probability level <P (+/3o)) are respectively solutions of the following 
optimization problems: 

q- = minF- 1 (*(-Ai)|0). i = l,-..,n (26) 
q+ =maxF x 1 (<P(+p )\9), i = l,...,n (27) 

where {F" 1 , i = l, ... ,n} are the quantile functions of the margins. If the domain @ e is rectangular and if one is 
able to derive an analytical expression for the quantile functions of the margins and their derivatives with respect to 
the parameters, then these optimization problems might be solved analytically. However assuming a more general 
setup where one has only numerical definition of these quantities, these problems can be efficiently solved by 
means of a simple gradient-based algorithm due to the convenient properties of the quantile functions - namely, 
the monotony with respect to the location and shape parameters. Finally, the sought hyperrectangle can be easily 
defined by means of the following indicator function: 

, . |1 if q.T < Vs < q,t, i = l, ... ,n 
Wv)={„ otherwise * (28) 

Again, the normalizing constant of this PDF could be easily derived (hyperrectangle volume) though it is not 
required by the refinement procedure proposed in Section 3. 



4.2 The adaptive surrogate-based nested RBDO algorithm 

4.2.1 Optimization 

The kriging surrogate together with its adaptive refinement procedure is finally plugged into a double-loop RBDO 
algorithm. The outer optimization loop is performed by means of the Polak-He optimization algorithm (Polak, 
1997). Provided an initial design, this algorithm proceeds iteratively in two steps: (i) the direction of optimization 
is determined solving a quasi-SQP sub-optimization problem and (ii) the step size is approximated by the Goldstein- 
Armijo approximate line-search rule. 

4.2.2 Surrogate-based reliability and reliability sensitivity analyses 

The nested reliability and reliability sensitivity analyses are performed with the subset simulation variance reduction 
technique (Au and Beck, 2001) onto the kriging surrogates. The subset simulation technique for reliability sensi- 
tivity analysis is detailed in Song et al. (2009). Briefly, it takes advantage of the definition of the failure probability 
given in Eq. (2) . Indeed, pointing out that the limit-state equation does not explicitly depend on the design vari- 
ables 9, the differentiation of the failure probability only requires the differentiation of the joint PDF which can be 
derived analytically when the probabilistic model is defined in terms of margin distributions and copulas (Lebrun 
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and Dutfoy, 2009). The trick is inspired from importance sampling and proceeds as follows: 



86 



— fxix,0)dx 

86 >gW<o 



dfx(x.O) 



g(x)<0 



86 



dx 



36 



= E 



3fx(x,e) 
(X) 36 



f x (x, 0)dx 



'fx 



f x (x,e) 

The latter quantity is known to have an unbiased consistent estimator which reads as follows: 



8P 



86 N 



N d fx {x«\e) 



k=l 



fx[x«\9) 



(29) 
(30) 

(31) 

(32) 



(33) 



where the sample {x 



(i) 



x^} is the same as the one used for the estimation of the failure probability P f . In 



sp, 



other words, the estimation of -jj does not require any additional simulation runs: it simply consists in a post- 
processing of the samples generated for reliability estimation. The concept can be easily extended to the subset 
simulation technique - see Song et al. (2009) for the details. 



4.2.3 Implementation 



The overall methodology was implemented within the FERUM v4.0 toolbox (Bourinet et al, 2009). It makes use 
of the Matlab toolboxes functions quadprog (for the SQP sub-optimization problem) and slicesample (to gen- 
erate samples from the refinement PDF t €). We provide a summarized pseudo-code of the proposed strategy in 
Algorithm 2. 



Algorithm 2 Adaptive surrogate-based nested RBDO strategy 



1: (o) , e L , B u ,j:=0 
2: for i = 1 to n do 

3: 9w ~ „™ n .„ p x,H x i> 9 ) 



L <e<e 



4: q+:= max F" 1 (x f , 0) 

1 e L <e<e u 

5: end for 

6: w := v ^ 1 %j/3o (v) 
7: k := 4>- 1 (97.5%), s p := 10" 1 
8: Refine := true, Optimize := true 
9: while Refine and Optimize do 

10 

11 

12 

13 

14 
15 
16 
17 
18 
19 
20 
21 

22 
23 
24 



25 
26 



while Refine do 

G := RefineKrigingModel(w, k, ep) <-» Use Algorithm 1 
end while 

F°:={x : ix 5 {x)<0} 

/3 0) , V fl /3 0) := ReliabilityAnalysis(F °, O) ) 
c t0 :=C (0O)) ; VflC 0) := v ,c(9(i)) 

/O) :=/(0 a) );Ve/ o) . = Vef(e m ) 

d.V> := SolveQuasiSQP(c«, f { '\ /8«, V e c<-», V e f { '\ V e P^) 
s 0) ■- GoldsteinArmijoStepSizeRule(c, /, /3, d 0) ) 
(j+i). = e O) +s O) d O) 

for£:=-l, 0,-1-1 do 

F ' -={x : n d (x) + ika s (x)<0} 

P 1 := ReliabilityAnalysis(F ; , 6 {i+1) ) 
end for 

Refine := max - /3 °, /3 - /3" 1 ) > e p 
Optimize := \\e { ' +1) - > s e or ||c 0+1) - c 0) || > s c or 3 i \ f t > or /? < / 



end while 
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The first step of the algorithm consists in finding the hyperrectangular region that bounds the confidence region 
of the augmented probability density function according to Section 4.1. Once this is done, one may define the 
uniform weighting density w = ly Vi p in Eq. (28) and use it within the adaptive population-based refinement 
procedure detailed in Section 3 and Algorithm 1. Kriging models are built for each performance function, and they 
are refined until they meet the selected accuracy regarding reliability estimation. It is worth noting that the kriging 
surrogates are built in the augmented reliability space spanned by V, but used in the current space of random 
variables spanned by X \ = As soon as they are accurate enough we perform surrogate-based reliability 
and reliability sensitivity analysis in order to propose an improved design. A quasi-SQP algorithm is then used in 
order to determine the best improvement direction; and the Goldstein-Armijo approximate line-search rule is used 
to find the best step size along that direction. The current design is improved and the kriging model accuracy for 
reliability estimation is being checked at the improved design. The convergence is obtained if the optimization has 
converged (using the regular criteria in gradient-based deterministic optimization) and if the kriging models allow 
a sufficiently accurate reliability estimation according to the proposed error measure. 



5 Applications 

In this section, the proposed adaptive nested surrogate-based RBDO strategy is applied to some examples from 
the literature for performance comparison purposes. All the kriging surrogates are sequentially refined in order to 
achieve an empirical error measure Sp < 1CT 1 on the estimation of the reliability indices. 



5.1 Elastic buckling of a straight column - An analytical reference 

In essence, the purpose of this first basic example is to validate the proposed algorithm with respect to a reference 
analytical solution. 

Let us consider a long simply-supported rectangular column with section bxh subjected to a constant service 
axial load F ser . Provided h < b and its constitutive material is characterized by a linear elastic behavior through its 
Young's modulus E, its critical buckling load is given by the Euler formula: 

n 2 Ebh 3 

F* = -^- (34) 
This allows one to formulate the performance function which will be involved in the probabilistic constraint as: 

n 2 Ebh 3 

g (£, b, h, L) = 12£2 - F ser . (35) 

The probabilistic model consists in the 3 independent random variables given in Table 1. 

Applying the performance function as a design rule in a fully deterministic fashion (i.e. using the means of the 
random variates) allows to determine the service load F ser so that the initial deterministic design /x^ = /jl^ = 
200 mm satisfies the limit-state equation g(x) = 0: 

*« = —$- = ^ • (36) 

The reliability-based design problem consists in finding the optimal means fi b and fj, h of the random width b 
and height h. The optimal design is the one that minimizes the average cross section area which is approximated 
as follows: 

c{0) = n b n h . (37) 
It should also satisfy the following deterministic constraint: 

f{8) = i x h - i x i <Q (38) 



Variable 


Distribution 


Mean 


C.o.V 


E (MPa) 


Lognormal 


10000 


15% 


b (mm) 


Lognormal 


Ms 


5% 


h (mm) 


Lognormal 


P-h 


5% 


L (mm) 


Deterministic 


3 000 





Table 1 Elastic buckling of a straight column - Probabilistic model 
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in order to ensure that the Euler formula is applicable, as well as the following safety probabilistic constraint: 

-#- 1 (P(g(X)<0)) >p (39) 
where /3 = 3 is the generalized target reliability index. 

5.1.2 Analytical solution 

Due to the use of lognormal random variates in the probabilistic model together with the simple performance func- 
tion (multiplications and divisions), the problem can be handled analytically. Indeed, the isoprobabilistic transform 
of the limit-state surface equation in terms of standard normal random variates is straightforward and leads to a 
linear equation. And it finally turns out after basic algebra that the Hasofer-Lind reliability index (associated with 
the exact failure probability) reads: 

1 °gfT^l +A fi + A b + 3A "-2A i 

/^hl = v ; '{ = = 2 (40) 

where X. and denote the parameters of the lognormal random variates. 

The optimal solution of the RBDO problem is then simply derived by saturating the two constraints in log-scale 
(i.e. with respect to X b and A h ) and this leads to the square cross section with parameters: 

K =K = \ (ft 4? E + e b + %\ + 4£f - log (^rj -A £ + 2A^. (41) 
5.1.2 Numerical solution 

The proposed numerical strategy is applied in order to solve the RBDO problem numerically. The refinement proce- 
dure of the limit-state surface is initialized with an initial DOE of K = 10 points and K = 10 points are sequentially 
added to the DOE if it is not accurate enough for reliability estimation. 




02468 02468 02468 

Iterations Iterations Iterations 



(a) Run #1 - Starting from the optimal deterministic design 




200 1 ' 1 5 1 1 1 1 ' ' 

024602460246 
Iterations Iterations Iterations 



(b) Run #2 - Starting from a conservative design 



Fig. 3 Elastic buckling of a straight column - Convergence 
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The convergence of the algorithm is depicted in Figure 3 for two runs starting from different initial designs. Run 
#1 is initiated with the optimal deterministic design [i b = fji h = b* = h* = 200 mm whereas Run #2 is initiated with 
an oversized design fx b = fx h = 300 mm. Convergence is achieved as all the constraints (deterministic and reliability- 
based) are satisfied and both the cost and design variables have reached a stable value. The algorithm converges 
to the exact solution derived in the previous subsection which is the square section with width ix* b = juf m 231 mm 
- the approximation of the exact solution is only due to the numerical error. Note that the reliability-based optimal 
design is 15% higher than the optimal deterministic design for the chosen reliability level (/3 = 3). 

This optimum is reached using only 20 evaluations of the performance function thanks to the kriging surrogate. 
The DOE used for this purpose is enriched only once and it is then accurate enough for all the design configurations 
including the optimal design. Running the same RBDO algorithm without using the kriging surrogates (i.e. using 
subset simulation onto the real performance function for the nested reliability and reliability sensitivity analyses) 
requires about 4 x 10 6 evaluations of the performance function for the same number of iterations of the optimizer 
and converges to the same optimal design. 

5.2 Short column case 

This simple mechanical example is extensively used in the RBDO literature as a benchmark for numerical methods 
in RBDO. In this paper, we use the results from the article by Royset et al. (2001) as reference. It consists in a short 
column with rectangular cross section b x h. It is subjected to an axial load F and two bending moments M 1 and 
M 2 whose axes are defined with respect to the two axes of inertia of the cross section. Such a load combination 
is referred to as oblique bending due to the rotation of the neutral axis. Under the additional assumption that the 
material is elastic perfectly plastic, the performance function describing the ultimate serviceability of the column 
with respect to its yield stress f y reads as follows: 

t \ 4Mj 4M 2 f F V 

G^M^.M)-!- — - — -^ — J . (42) 

The stochastic model originally involves three independent random variables whose distributions are given in 
Table 2. Note that in the original paper, the design variables [i b and /j, h are considered as deterministic. Since the 
present approach only deals with design variables that defines the joint PDF of the random vector X, they are 
considered here as Gaussian with a small coefficient of variation and the optimization is performed with respect to 
their mean [i h and [i h . The objective function is formulated as follows: 

c(Mb. Phi = c iHb, Mh) +Pf(Hb, Hh)c f (Hb, Vh) = (1 + 100 P f ) (43) 

where the product Pf x Cf is the expected failure cost which is chosen as proportional to the construction cost 
c . The search for the optimal design is limited to the designs that satisfy the following geometrical constraints: 
1/2 < ix b /ix h < 2 with 100 < fj, b , fj, h < 1000, and the minimum reliability is chosen as p = 3. 



Variable 


Distribution 


Mean 


C.o.V 


Mj (N.mm) 


Lognormal 


250 x 10 6 


30% 


M 2 (N.mm) 


Lognormal 


125 x 10 6 


30% 


P (N) 


Lognormal 


2.5 x 10 6 


20% 


f y (MPa) 


Lognormal 


40 


10% 


b (mm) 


Gaussian 


f*b 


1%* 


h (mm) 


Gaussian 




1%* 



* Additional uncertainty introduced to match the RBDO problem formulation in Eq. (1) 



Table 2 Short column case from Royset et al. (2001) - Probabilistic model 



The results are given in Table 3. In this table, /3 HL denotes the Hasofer-Lind reliability index (FORM-based), 
and /3 sim denotes the generalized reliability index estimated by subset simulation with a coefficient of variation less 
than 5%. The deterministic design optimization (DDO) was performed using the mean values of all the variables 
in Table 2 without considering uncertainty and thus leads to a 50% failure probability. Note that the corresponding 
optimal cost does not account for the expected failure cost. The other lines of Table 3 shows the results of the 
RBDO problem. The first row gives the reference results from Royset et al. (2001). The number of performance 
function calls was not given in the original paper. However it may be estimated to 4 x 10 6 given the methodology 
the authors used and assuming they targeted a 5% coefficient of variation on the failure probability in their Monte- 
Carlo simulation. The second row provides the results from a FORM-based nested RBDO algorithm (RIA). This 
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Method 


Opt. design (mm) 


Cost (mm 2 ) 


# of func. calls 


Reliability 


DDO 


b = 258 h = 500 


1.29 x 10 b 


50 


ftim « 0.01 


RBDO OS > 3) 










Reference (DSA) 


b = 372 h = 559 


2.15 x 10 5 




ftim 3-38 


FORM-based (RIA) 


b = 399 7l=513 


2.12 x 10 5 


9472 


ftn. = 3.38 


Present w/o kriging 


b = 369 h = 560 


2.16 x 10 5 


19 x 10 6 


ftim 3-35 


Present w/ kriging 


b = 379 h = 547 


2.17 x 10 5 


140 


ftim 3-32 



Table 3 Short column case from Royset et al. (2001) - Comparative results 



latter approach seems to lead to a slightly better design though it is due to the first-order reliability assumptions 
that are not conservative in this case. Indeed, subset simulation leads to a little lower generalized reliability index 
ftim 3.19 (with a 5% C.o.Y) which in turns slightly increases the failure-dependent objective function to c « 
2.20 x 10 5 mm 2 . This example shows that FORM-based approaches can mistakenly lead to non conservative optimal 
designs - without any self-quantification of the possible degree of non conservatism. The third row gives the results 
obtained by the same nested RBDO algorithm, using however the subset simulation technique as the reliability (and 
reliability sensitivity) estimator. Finally, the last row gives the results obtained when using kriging as a surrogate for 
the limit-state surface. The kriging model used for this application used a constant regression model and a squared 
exponential autocovariance model. It was initialized with a 50-point DOE and sequentially refined with K = 10 
points per refinement iteration. 




Fig. 4 Short column case - Convergence 



Another interesting fact about this example is that the reliability constraint is not saturated at the optimum: the 
algorithm converges at a higher reliability level as illustrated in Figure 4. This is due to the specific formulation of 
the cost function in Eq. (43) that accounts for a failure cost that is indexed onto the failure probability. Indeed, the 
cost function behaves itself as a constraint and the optimal reliability level is formulated in terms of an acceptable 
risk (probability of occurrence times consequence) instead of an acceptable reliability index /3 . 



5.3 Bracket structure case 

This mechanical example is originally taken from Chateauneuf and Aoues (2008) . It consists in the study of the 
failure modes of the bracket structure pictured in Figure 5. The bracket structure is loaded by its own weight due 
to gravity and by an additional load at the right tip. The two failure modes under considerations are: 

- the maximum bending in the horizontal beam (CD, at point B) should not exceed the yield strength of the 
constitutive material, so that the first performance function reads as follows: 

gi (w C d, t, L, P, p, fy) =f y -cT B (44) 
where the maximum bending stress reads: 



6M B PL pgw CD tL 



2 



with: Mg = — — — h " , (45) 



''CD 



t 2 ~'~ D 3 ' 18 
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Fig. 5 Bracket structure case - Mechanical model 



- the maximum axial load in the inclined member (AB) should not exceed the Euler critical buckling load (ne- 
glecting its own weight), so that the second performance function reads as follows: 

g 2 (waj,, w cd , t, L, P, p, / y ) = F buckling - Fab (46) 

where the critical Euler buckling load F buckling is defined as: 

n 2 EI n 2 E t w 3 ^ 9 sin 2 6 

^buckling — ^2 — 48L 2 ' ^ 

and the resultant of axial forces in member AB reads (neglecting its own weight) : 

1 (3P 3pgw CD tL } 

F ^ = ^e{T + 4 J- (48) 

The probabilistic model for this example is the collection of independent random variables given in Table 4. Note 
that the coefficient of variation of the random design variables is kept constant along the optimization as in the 
original paper. 



Variable 


Distribution 


Mean 


C.o.V 


P 


(kN) 


Gumbel 


100 


15% 


E 


(GPa) 


Gumbel 


200 


8% 


fy 


(MPa) 


Lognormal 


225 


8% 


P 


(kg/m 3 ) 


Weibull 


7860 


10% 


L 


(m) 


Gaussian 


5 


5% 


Wab 


(mm) 


Gaussian 




5% 


wqd 


(mm) 


Gaussian 


Mwcd 


5% 


t 


(mm) 


Gaussian 


n t 


5% 



Table 4 Bracket structure case - Probabilistic model 



The RBDO problem consists in finding the rectangular cross sections of the two structural members that mini- 
mize the expected overall structural weight which is approximated as follows: 

c (m Wab > Mw CD ,t)=ix pi x ti x L ( ^y- ix Wm + MwcD J (49) 

while satisfying a minimum reliability requirement equal to /3 = 2 with respect to the two limit-states in Eq. (44) 
and Eq. (46) considered independently. The search for the optimal design is bounded to the following hyperrectan- 
gle: 50 < fJL WAs , [i WcD , fit ^ 3 °0 (in mm). 

The comparative results are provided in Table 5 considering the results from Chateauneuf and Aoues (2008) as 
reference. 
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Method Opt. design (mm) Cost (kg) # of func. calls Reliability 



DDOw/PSF* w CD = 202 2 632 40 ' siml 

t = 269 ,sim2 

RBDO (ft > 2, i = 1, 2) 

Wab = 61 
SORA* w CD = 157 

t = 209 



RIA* 


W C D 
t 


= 61 
= 157 
= 209 


1675 


2 340 


ftiml 
ftim2 


ss 1.96 
sb 2.01 


Present w/o kriging 


W C D 
t 


= 58 
= 119 
= 241 


1550 


5 x 10 6 
5 x 10 6 


ftiml 
ftim2 


ss 2.00 
ss 2.01 


Present w/ kriging 


wab 
w CD 
t 


= 59 
= 135 
= 226 


1610 


100 
150 


ftiml 
ftim2 


ss 2.01 
ss 2.03 



* As computed by Chateauneuf and Aoues (2008). 



Table 5 Bracket structure case - Comparative results 



1675 



1340 



ftiml 
ftim2 



1.96 
2.01 



The first row gives the deterministic optimal design that was obtained by Chateauneuf and Aoues (2008) using 
partial safety factors (PSF). It can be seen from the reliability indices that these PSF provide a significant safety 
level. However, one may rather want to find an even lighter design allowing for a lower safety level /3 = 2. To 
do so, the RBDO formulation of the problem is solved. Chateauneuf and Aoues (2008) used the SORA technique 
which is a decoupled FORM-based approach. The reliability indices at the optimal design were checked using 
the subset simulation technique (targeting a coefficient of variation less than 5%) and revealed that the FORM- 
based approach slightly underestimates the first optimal reliability index in this case. The RIA technique which is a 
standard double-loop FORM-based approach provides the same solution but it is less efficient. 

Implementing the proposed approach without plugging the kriging surrogates converges to similar results but 
clearly confirms that direct simulation-based approaches are not tractable for RBDO. Replacing the performance 
function by their kriging counterparts allows to save a significant number of simulations (10 2 opposed to 10 6 ) and 
in addition, to provide an error measure on the reliability estimation as opposed to the FORM-based approaches. 
However, one may note the disparities between the proposed designs. First, the disparity between FORM-based 
methods and the presently proposed strategies is explained by the conservatism of the FORM assumptions in this 
case. Then, the disparity between the two present approaches is certainly due to the flatness of the sub-optimization 
problem and the stochastic nature of the simulation-based reliability estimation. The convergence of the algorithm 
is depicted in Figure 6. 




Fig. 6 Bracket structure case - Convergence 
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6 Conclusion 

The aim of the present paper was to develop a strategy for solving reliability-based design optimization (RBDO) 
problems that is applicable to engineering problems involving time-consuming performance models. Starting with 
the premise that simulation-based approaches are not affordable when the performance function involves the out- 
put of an expensive-to-evaluate computer model, and that the MPFP-based approaches do not allow to quantify the 
error on the estimation of the failure probability an approach based on kriging and subset simulation is explored. 

The strategy has been tested on a set of examples from the RBDO literature and proved to be competitive with 
respect to its FORM-based counterparts. Indeed, convergence is achieved with only a few dozen evaluations of the 
real performance functions. In contrast with the FORM-based approaches, the proposed error measure allows one 
to quantify and sequentially minimize the surrogate error onto the final quantity of interest: the optimal failure 
probability. It is important to note that the numerical efficiency of the proposed strategy mainly relies upon the 
properties of the space where the kriging surrogates are built: the so-called augmented reliability space. This space 
is obtained by considering that the design variables in the RBDO problem simply augments the uncertainty in the 
random vector involved in the reliability problem. Building the surrogates in such a space allows one to reuse 
them from one RBDO iteration to the other and thus saves a large number of performance functions evaluations. 
It is also worth noting that the original refinement strategy proposed in Section 3 makes it possible to add several 
observations in the design of experiments at the same time, and thus to benefit from the availability of a distributed 
computing platform to speed up convergence. 

However, as already mentioned in the literature, it was observed that the number of experiments increases with 
the number of variables involved in the performance functions and that the kriging strategy loses numerical effi- 
ciency when the DOE contains more than a few thousands experiments - although such an amount of information 
is not even available in real-world engineering cases. This latter point requires further investigation. A problem 
involving a nonlinear-finite-element-based performance function and 10 variables is currently investigated and will 
be published in a forthcoming paper. 
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